# -*- coding: utf-8 -*-
"""
Created on Thu Jul 28 08:18:36 2022

@author: Derek Rodriguez Pacheco, PhD Candidate
IUSS Pavia
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import subprocess
import matplotlib as mpl
from Spectra.SpectrumA import spectrumAD 

mpl.rc('font',family = 'Times New Roman')

m= 970                           ## Mass in N*s2/m or kg
g= 9.81                          ## Acceleration of gravity in m/s2

adjfactors= np.array(pd.read_csv('Adjustment_Factors_Displacement_CM.txt', header= None, delim_whitespace= True)[0])
indexesb= np.array(pd.read_csv('Indexes_Beginning.txt', header= None, delim_whitespace= True)[0])
indexesf= np.array(pd.read_csv('Indexes_End.txt', header= None, delim_whitespace= True)[0])

## Model's parameters

k1= 0.5
k2= 0.1
sigAct= 3.4
betaf= 1
epsSlip= 0
epsBear= 0
rBear= 1

afs= 0.8
afd= 0.8


# s1p= 3*afs
s1p= 2.5
# e1p= 2*afd
e1p= 6
# s2p= 9*afs
s2p= 4
# e2p= 10*afd
e2p= 25
# s3p= 35*afs
s3p= 15
# e3p= 55*afd
e3p= 50
pinchX= 0.5
pinchY= 0.01
damage1= 0
damage2= 0
beta= 0

s1n= s1p*45
e1n= np.copy(e1p)
s2n= s2p*45
e2n= np.copy(e2p)
s3n= s3p*45
e3n= np.copy(e3p)

width= 800
height= 1076

testnumber= [13,14,15,16,17,18,19,20,24,25,26,27,28,40,41]

##

datacm= {}
# acceltcm01= {}
# acceltcm02= {}
accelcm= {}
timecm= {}
timecmc= {}

data= {}
datad= {}
accelda= {}
# accelda1= {}
# accelro= {}
accelbase= {}
dispda= {}
dispdacm= {}
timearray= {}
energyarray= {}
energyarraycm= {}
maxaccel= []
minaccel= []
maxaccelf= []
matchindexmin= []
matchindexmax= []
matchdispmin= []
matchdispmax= []
timetotal= []
acceltotal= []
dispnet= {}
dispnetcm= {}
baserot = {}
# timestart= []

for i in range(15):
    datacm['%d' % (i+1)]= np.array(pd.read_csv('Data_Center_Mass\\accelAndDispl_15Hz\\%d_AccelAndDispl_15Hz.txt' % testnumber[i], header= None, delim_whitespace= True))
    accelcm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,4])[indexesb[i]:indexesf[i]+1]
    timecm['%d' % (i+1)]= np.array(datacm['%d' % (i+1)][:,0])[indexesb[i]:indexesf[i]+1]
    timecmc['%d' % (i+1)]= timecm['%d' % (i+1)]-timecm['%d' % (i+1)][0]

data113= np.array(pd.read_csv('13_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
data['1']= np.array(pd.read_csv('13_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['2']= np.array(pd.read_csv('14_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['3']= np.array(pd.read_csv('15_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['4']= np.array(pd.read_csv('16_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['5']= np.array(pd.read_csv('17_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['6']= np.array(pd.read_csv('18_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['7']= np.array(pd.read_csv('19_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['8']= np.array(pd.read_csv('20_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['9']= np.array(pd.read_csv('24_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['10']= np.array(pd.read_csv('25_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['11']= np.array(pd.read_csv('26_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['12']= np.array(pd.read_csv('27_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['13']= np.array(pd.read_csv('28_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['14']= np.array(pd.read_csv('40_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))
data['15']= np.array(pd.read_csv('41_time_baseAcc_displGross_displNet_topACC_filterLP_15Hz.txt', header= None, delim_whitespace= True))

datad['1']= np.array(pd.read_csv('13_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['2']= np.array(pd.read_csv('14_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['3']= np.array(pd.read_csv('15_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['4']= np.array(pd.read_csv('16_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['5']= np.array(pd.read_csv('17_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['6']= np.array(pd.read_csv('18_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['7']= np.array(pd.read_csv('19_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['8']= np.array(pd.read_csv('20_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['9']= np.array(pd.read_csv('24_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['10']= np.array(pd.read_csv('25_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['11']= np.array(pd.read_csv('26_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['12']= np.array(pd.read_csv('27_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['13']= np.array(pd.read_csv('28_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['14']= np.array(pd.read_csv('40_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))
datad['15']= np.array(pd.read_csv('41_t_rockRot_rockDispl_grossDispl_netDispl_netENE.txt', header= None, delim_whitespace= True))

for i in range(15):
    accelda['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,4])
    accelbase['%d' % (i+1)]= np.array(data['%d' % (i+1)][:,1])
    dispda['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,3])
    dispnet['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,4])
    dispdacm['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,3])*adjfactors[i]
    dispnetcm['%d' % (i+1)]= np.array(datad['%d' % (i+1)][:,4])*adjfactors[i]
    timearray['%d' % (i+1)]= data['%d' % (i+1)][:,0]-data['%d' % (i+1)][0,0]
    energyarray['%d' % (i+1)]= [np.trapz(accelda['%d' % (i+1)][:j]*m,dispda['%d' % (i+1)][:j]/1000) for j in range(len(accelda['%d' % (i+1)]))]
    energyarraycm['%d' % (i+1)]= [np.trapz(accelcm['%d' % (i+1)][:j]*m,dispdacm['%d' % (i+1)][:j]/1000) for j in range(len(accelcm['%d' % (i+1)]))]
    baserot['%d' % (i+1)] = datad['%d' % (i+1)][:,1] 
cc= 0

for i in range(15):
    for j in range(len(timearray['%d' % (i+1)])):
        if i==0:
            timetotal.append(timearray['%d' % (i+1)][j])
        else:
            timetotal.append(cc+timearray['%d' % (i+1)][j])
        acceltotal.append(accelbase['%d' % (i+1)][j])
    cc= timetotal[-1]
    if i!=14:
        for j in range(int(round(15/timearray['%d' % (i+1)][1]))):
            if i==0:
                timetotal.append(timearray['%d' % (i+1)][-1]+(j+1)*timearray['%d' % (i+1)][1])
            else:
                timetotal.append(cc+(j+1)*timearray['%d' % (i+1)][1])
            acceltotal.append(0)
        cc= timetotal[-1]

#### SDOF Analysis

# for i in range(15):
np.savetxt('TH/Time_Total.txt', timetotal)
np.savetxt('TH/Acceleration_Total.txt', acceltotal)

recordsa= 'TH/Acceleration_Total.txt'
recordst= 'TH/Time_Total.txt'

dts= 0.0010
mass= m/1e6

dur= timetotal[-1]

fh= open('FMdata.tcl', 'w')

fh.write(
    'set FMfilet "%s";\nset FMfilea "%s";\nset dtg %f;\nset durs %f;'
    '\nset M %f;'
    '\nset s1p %f;\nset e1p %f;\nset s2p %f;\nset e2p %f;\nset s3p %f;'
    '\nset e3p %f;\nset pinchX %f;\nset pinchY %f;\nset damage1 %f;'
    '\nset damage2 %f;\nset beta %f;' 
    '\nset s1n %f;\nset e1n %f;\nset s2n %f;\nset e2n %f;\nset s3n %f;\nset e3n %f;'
    '\nset width %d;\nset height %d;'
    '\nset k1 %f;\nset k2 %f;\nset sigAct %f;\nset betaf %f;\nset epsSlip %f;\nset epsBear %f;\nset rBear %f;' % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height, k1, k2, sigAct, betaf, epsSlip, epsBear, rBear))
    # % (recordst, recordsa, dts, dur, mass, s1p, e1p, s2p, e2p, s3p, e3p, pinchX, pinchY, damage1, damage2, beta, s1n, e1n, s2n, e2n, s3n, e3n, width, height))
fh.close()
subprocess.call('OpenSees MDOF_CalibrationMH_Grav.tcl', shell=True)

####

sdofd= {}
sdofa= {}
sdoft= {}
sdoff= {}
sdofe= {}
sdofr= {}

indexes= [0,43014,86028,129057,171918,214781,257642,300501,343403,386241,427407,468579,510116,551919,594286,622165]
timei= [0,43.014,86.028,129.057,171.918,214.781,257.642,300.501,343.403,386.241,427.407,468.579,510.116,551.919,594.286]

for i in range(15):
    sdofd['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/DispT.out', header= None, delim_whitespace= True)[1][indexes[i]:indexes[i+1]])
    sdofa['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[1][indexes[i]:indexes[i+1]])
    sdoft['%d' % (i+1)]= np.array(pd.read_csv('ResultsAtt/AccelT.out', header= None, delim_whitespace= True)[0][indexes[i]:indexes[i+1]])-timei[i]
    sdoff['%d' % (i+1)]= sdofa['%d' % (i+1)]*0.001*m*-1
    sdofe['%d' % (i+1)]= [np.trapz(sdoff['%d' % (i+1)][:j],sdofd['%d' % (i+1)][:j]*0.001) for j in range(len(sdoff['%d' % (i+1)]))]
    upR = np.array(pd.read_csv('ResultsAtt/UpliftT.out', header= None, delim_whitespace= True)[2][indexes[i]:indexes[i+1]])
    upL = np.array(pd.read_csv('ResultsAtt/UpliftT.out', header= None, delim_whitespace= True)[1][indexes[i]:indexes[i+1]])
    sdofr['%d' % (i+1)] = (upL-upR)/width


## Plotting

testnum2= [2,3,4,5,6,7,8,9,11,12,13,14,15,17,18]
rtzero= [28013,28013,28028,27860,27860,27860,27860,27900,27840,27170,26170,26540,26800,27370,27878]
intensity= [0.17,0.17,0.19,0.20,0.20,0.19,0.20,0.40,0.29,0.29,0.34,0.39,0.46,0.51,0.52]
location1= [2500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500,3500]
location2= [4000,11000,11000,11000,11000,11000,11000,11000,11000,11000,11000,11000,11000,11000,11000]
location3= [-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58,-58]
location4= [5.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0,7.0]
subplotnum = ['a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p']
nr= 3
nc= 5

fig1, ax1 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        ax1[i,j].plot(timecmc['%d' % count], energyarraycm['%d' % count], color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax1[i,j].plot(sdoft['%d' % count][:rtzero[count-1]], sdofe['%d' % count][:rtzero[count-1]], color= '#FF8787', lw= 1.5, label= 'Numerical model')
                
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax1[i,j].set_xlabel('Time (s)')
            # ax1[i,j].set_xticks([0,10,20,30])
            # ax1[i,j].set_xlim([0,30])
        if count==1 or count==6 or count==11:
            ax1[i,j].set_ylabel('Absorbed Energy (J)')
            # ax1[i,j].set_yticks([0,800,1600,2400,3200])
            # ax1[i,j].set_ylim([0,3200])
        ax1[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, location1[count-1]), 
                      xycoords = 'data',
                      xytext = (0, -5), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax1[i,j].annotate('%s)' % subplotnum[count-1], xy=(27,3700),fontsize=14)
        ax1[i,j].set_xticks([0,10,20,30])
        ax1[i,j].set_yticks([0,1000,2000,3000,4000])
        ax1[i,j].set_xlim([0,30])
        ax1[i,j].set_ylim([0,4000])
        if count!=1 and count!=6 and count!=11:
            ax1[i,j].set_yticklabels([])
        if count!=11 and count!=12 and count!=13 and count!=14 and count!=15:
            ax1[i,j].set_xticklabels([])
        ax1[i,j].set_facecolor('#DEE1E6')
        ax1[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax1[i,j].legend(loc=2)
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig1.tight_layout()        
#fig1.savefig('Absorbed_Energy.png', bbox_inches= 'tight',dpi=300)
fig1.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure11.tiff', bbox_inches= 'tight',dpi=300)


fig2, ax2 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        
        paccE = np.argmax(abs(accelcm['%d' % count]))
        paccN = np.argmax(abs(sdoff['%d' % count]))
        
        ax2[i,j].plot(dispdacm['%d' % count], accelcm['%d' % count]*m, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax2[i,j].plot(sdofd['%d' % count][:rtzero[count-1]], sdoff['%d' % count][:rtzero[count-1]], color= '#FF8787', lw= 1.5, label= 'Numerical model')
                
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax2[i,j].set_xlabel('Displacement (mm)')
        if count==1 or count==6 or count==11:
            ax2[i,j].set_ylabel('Force (N)')
        ax2[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (-58, location2[count-1]), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        if(count==1):
            y = -86
        else:
            y = -120
        ax2[i,j].annotate('Relative error of peak = '+str(round(((abs(sdoff['%d' % count][paccN])-abs(m*accelcm['%d' % count][paccE]))/abs(m*accelcm['%d' % count][paccE]))*100,1))+'%',
                      xy = (-58, location2[count-1]-1000), 
                      xycoords = 'data',
                      xytext = (+0, y), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax2[i,j].set_xticks([-60,-30,0,30,60])
        ax2[i,j].set_yticks([-15000,-10000,-5000,0,5000,10000,15000])
        ax2[i,j].set_xlim([-60,60])
        ax2[i,j].set_ylim([-15000,15000])
        ax2[i,j].annotate('%s)' % subplotnum[count-1], xy=(50,12700),fontsize=14)
        if count!=1 and count!=6 and count!=11:
            ax2[i,j].set_yticklabels([])
        if count!=11 and count!=12 and count!=13 and count!=14 and count!=15:
            ax2[i,j].set_xticklabels([])
        ax2[i,j].set_facecolor('#DEE1E6')
        ax2[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax2[i,j].legend(loc=2)
        
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig2.tight_layout()
fig2.savefig('Hysteretic_Response.png', bbox_inches= 'tight',dpi=300)        
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig2.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure13.tiff', bbox_inches= 'tight',dpi=300)

fig3, ax3 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        pcrdE = np.argmax(abs(dispdacm['%d' % count]))
        pcrdN = np.argmax(abs(sdofd['%d' % count]))
        
        ax3[i,j].plot(timecmc['%d' % count], dispdacm['%d' % count], alpha= 1, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax3[i,j].plot(sdoft['%d' % count][:rtzero[count-1]], sdofd['%d' % count][:rtzero[count-1]], alpha= 0.75, color= '#FF8787', lw= 1.5, label= 'Numerical model')
        #ax3[i,j].plot(timecmc['%d' % count][pcrdE],dispdacm['%d' % count][pcrdE],'o',mec='limegreen',mfc='none')
        #ax3[i,j].plot(sdoft['%d' % count][pcrdN],sdofd['%d' % count][pcrdN],'o',mec='#FF8787',mfc='none')
                
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax3[i,j].set_xlabel('Time (s)')
        if count==1 or count==6 or count==11:
            ax3[i,j].set_ylabel('Displacement (mm)')
        ax3[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, location3[count-1]+6), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax3[i,j].annotate('Relative error of peak = '+str(round(((abs(sdofd['%d' % count][pcrdN])-abs(dispdacm['%d' % count][pcrdE]))/abs(dispdacm['%d' % count][pcrdE]))*100,1))+'%',
                      xy = (1, location3[count-1]), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax3[i,j].annotate('%s)' % subplotnum[count-1], xy=(27,50),fontsize=14)
        ax3[i,j].set_xticks([0,10,20,30])
        ax3[i,j].set_yticks([-60,-40,-20,0,20,40,60])
        ax3[i,j].set_xlim([0,30])
        ax3[i,j].set_ylim([-60,60])
        if count!=1 and count!=6 and count!=11:
            ax3[i,j].set_yticklabels([])
        if count!=11 and count!=12 and count!=13 and count!=14 and count!=15:
            ax3[i,j].set_xticklabels([])
        ax3[i,j].set_facecolor('#DEE1E6')
        ax3[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax3[i,j].legend(loc=2)
        
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig3.tight_layout()    
fig3.savefig('Relative_Displacement.png', bbox_inches= 'tight',dpi=300)    
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig3.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure15.tiff', bbox_inches= 'tight',dpi=300)



fig4, ax4 = plt.subplots(nr,nc, figsize= (13,8))

            
for i in range(nr):
    for j in range(nc):
        count= nc*i+(j+1)
        pbrotE = np.argmax(abs(baserot['%d' % count]))
        pbrotN = np.argmax(abs(sdofr['%d' % count]))
        
        ax4[i,j].plot(timearray['%d' % count], baserot['%d' % count], alpha= 1, color= 'limegreen', lw= 1.5, label= 'Experimental')
        ax4[i,j].plot(sdoft['%d' % count], -sdofr['%d' % count], alpha= 0.75, color= '#FF8787', lw= 1.5, label= 'Numerical model')
        ax4[i,j].plot((0,30),(0.0018,0.0018),ls='--',color='k',lw=0.5,label=r'$0.01 \alpha$')
        ax4[i,j].plot((0,30),(-0.0018,-0.0018),ls='--',color='k',lw=0.5)
        #ax4[i,j].plot(timearray['%d' % count][pbrotE],baserot['%d' % count][pbrotE],'o',mec='limegreen',mfc='none')
        #ax4[i,j].plot(sdoft['%d' % count][pbrotN],-sdofr['%d' % count][pbrotN],'o',mec='#FF8787',mfc='none')
        
        
        if count==11 or count==12 or count==13 or count==14 or count==15:
            ax4[i,j].set_xlabel('Time (s)')
        if count==1 or count==6 or count==11:
            ax4[i,j].set_ylabel('Base rotation (rad)')
        ax4[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[count-1],intensity[count-1]),
                      xy = (1, -0.06), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax4[i,j].annotate('Relative error of peak = '+str(round(((abs(sdofr['%d' % count][pbrotN])-abs(baserot['%d' % count][pbrotE]))/abs(baserot['%d' % count][pbrotE]))*100,1))+'%',
                      xy = (1, -0.07), 
                      xycoords = 'data',
                      xytext = (+0, +0), 
                      textcoords = 'offset points', 
                      fontsize = 10,
                      # weight= 'bold',
                      color = 'black',
                      )
        ax4[i,j].set_xticks([0,10,20,30])
        ax4[i,j].set_yticks([-0.08,-0.06,-0.04,-0.02,0,0.02,0.04,0.06,0.08])
        ax4[i,j].set_xlim([0,30])
        ax4[i,j].set_ylim([-0.08,0.08])
        ax4[i,j].set_facecolor('#DEE1E6')
        ax4[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        if count==1:
            ax4[i,j].legend(loc=2)
        ax4[i,j].annotate('%s)' % subplotnum[count-1], xy=(26,97),fontsize=14)
        # ax1[i,j].set_xticks(np.arange(0.05, 0.75+0.10, 0.10))
        # ax1[i,j].set_xlim(0.05, 0.75)            
                
        # ax1[i,j].set_title(r'$q_a= %d, T_a/T_1$= 0' % count)
        
fig4.tight_layout()    
#fig3.savefig('MDOF_Relative_Displacement.png', bbox_inches= 'tight',dpi=300)    
# fig1.savefig('Fragility_Seismic_Demand_Paper_Mod.png', dpi= 600, bbox_inches= 'tight')
fig4.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/FigureBR2.tiff', bbox_inches= 'tight',dpi=300)







'''
dispmaxE = np.zeros(len(testnumber))
dispmaxN = np.zeros(len(testnumber))

VmaxE = np.zeros(len(testnumber))
VmaxN = np.zeros(len(testnumber))


for i in range(len(testnumber)):
    dispmaxE[i] = max(abs(dispdacm['%d' % (i+1)]))
    dispmaxN[i] = max(abs(sdofd['%d' % (i+1)]))*1000
    VmaxE[i] = max(abs(accelcm['%d' % (i+1)]*m))
    VmaxN[i] = max(abs(sdoff['%d' % (i+1)]))
    print(str(dispmaxN[i]/dispmaxE[i])+',')
'''   

Ta = np.linspace(0,1,100)

SaE = {}
SdE = {}
SaN = {}
SdN = {}

for i in range(len(testnumber)):
    SaE['%d' % (i+1)],SdE['%d' % (i+1)] = spectrumAD(timecmc['%d' % (i+1)],accelcm['%d' % (i+1)],Ta)
    SaN['%d' % (i+1)],SdN['%d' % (i+1)] = spectrumAD(sdoft['%d' % (i+1)],sdofa['%d' % (i+1)]/1000,Ta)


fig5, ax5 = plt.subplots(nr,nc,figsize=(13,8))

for i in range(nr):
    for j in range(nc):
        ns= nc*i+(j+1)
        relES = np.mean(abs(SaN['%d' % ns][20:]-SaE['%d' % ns][20:])/SaE['%d' % ns][20:]*100)
        ax5[i,j].set_facecolor('#DEE1E6')
        ax5[i,j].grid(color= 'black', ls= '--', lw= 0.1)
        ax5[i,j].plot(Ta,SaE['%d' % ns],color='limegreen',lw=1.5,label='Experiment')
        ax5[i,j].plot(Ta, SaN['%d' % ns],color='#FF8787', lw= 1.5, label= 'Numerical model')
        ax5[i,j].annotate('Test Number: %d\nPTA= %.2fg' % (testnum2[ns-1],intensity[ns-1]),
                     xy = (0.02, location4[ns-1]), 
                     xycoords = 'data',
                     xytext = (+0, -5), 
                     textcoords = 'offset points', 
                     fontsize = 10,
                     # weight= 'bold',
                     color = 'black',
                     )
        ax5[i,j].annotate('Mean relative error = '+str(round(relES,1))+'%',
                     xy = (0.02, location4[ns-1]-0.85), 
                     xycoords = 'data',
                     xytext = (+0, +0), 
                     textcoords = 'offset points', 
                     fontsize = 10,
                     # weight= 'bold',
                     color = 'black',
                     )
        if ns == 1 or ns == 6 or ns==11:
            ax5[i,j].set_ylabel('Spectral Acceleration (g)')
        if ns == 11 or ns == 12 or ns==13 or ns==14 or ns==15:
            ax5[i,j].set_xlabel('Period (s)')
        if ns == 1:
            ax5[i,j].legend(loc=2)
        
        if ns==13:
            ax5[i,j].set_xticks([0,0.2,0.4,0.6,0.8,1])
            ax5[i,j].set_yticks([0,1,2,3,4,5,6,7,8])
            ax5[i,j].set_xlim([0,1])
            ax5[i,j].set_ylim([0,8])
        else:
            ax5[i,j].set_xticks([0,0.2,0.4,0.6,0.8,1])
            ax5[i,j].set_yticks([0,1,2,3,4,5,6,7,8])
            ax5[i,j].set_xlim([0,1])
            ax5[i,j].set_ylim([0,8])
        
        ax5[i,j].annotate('%s)' % subplotnum[ns-1], xy=(0.88,7.2),fontsize=14)


fig5.tight_layout()   

#plt.savefig('MDOF_InSpectraAcc.tiff',dpi=300)
fig5.savefig('C:/Users/rober/Documents/ROSE/CADS/ElectricalCabinetPaper/Figure17.tiff', bbox_inches= 'tight',dpi=300)





